# ifndef _HEART_H_
# define _HEART_H_

# include"Bspline.h"
# include<fstream>

// Calculate sqrt(3) and restore it.
const double sqr3 = sqrt(3);

// Function x = x(t)
class xfun : public Function
{
public:
    double operator() (const double & t) const
    {
        return (sqr3 * sin(2.0 * M_PI * t));
    }
};

// Function y = y(t)
class yfun : public Function
{
public:
    double operator() (const double & t) const
    {
        double k = 2.0 * M_PI * t;
        return (2.0 / 3.0 * (sqrt(fabs(sin(k)) * sqr3) + sqr3 * cos(k)));
    }
};

// The class can get interpolation curve by different numbers of knots.
class heart
{
private:
    int n;                  // Different kinds of interpolation numbers.
    Bspline Bx, By;         // Construct the function (0,1) to R2
    double length;          // The length of curve

    // Calculate norm-2
    double norm2(const struct point & s1, const struct point & s2)
    {
        return (sqrt((s1.x - s2.x) * (s1.x - s2.x) + (s1.y - s2.y) * (s1.y - s2.y)));
    };

public:
    heart() {};
    
    heart(const int & _n) : n{_n}
    {
        // Only can construct by even numbers and n >=8.
        if(n % 2 != 0)
        {
            std::cout << "The number of node need to be even." << std::endl;
            return;
        }
        if(n < 8)
        {
            std::cout << "The knots are not enough." << std::endl;
            return;
        }

        // Construct the knots
        std::vector<struct point> knot;
        int n2 = n / 2;
        double eps = 1.0 / n;
        xfun fx; yfun fy;
        for(int i = 0; i < n; i++)
        {
            if(i == 0)
            {
                struct point p0{0, 2.0 / sqr3};
                knot.push_back(p0);
            }
            else if(i == n2)
            {
                struct point pn2{0, -2.0 / sqr3};
                knot.push_back(pn2);
            }
            else if(i == 1 || i == n2 + 1)
            {
                struct point p1{fx(i * eps - eps * 2.0 / 3.0), fy(i * eps - eps * 2.0 / 3.0)};
                knot.push_back(p1);
            }
            else if(i == n-1 || i == n2 - 1)
            {
                struct point pn{fx(i * eps + eps * 2.0 / 3.0), fy(i * eps + eps * 2.0 / 3.0)};
                knot.push_back(pn);
            }
            else
            {
                struct point pi{fx(i * eps), fy(i * eps)};
                knot.push_back(pi);
            }
        }

        // Comput cumulative chordal lengths.
        std::vector<double> tpoint;
        tpoint.push_back(0);
        for (int i = 1; i < n; i++)
        {
            tpoint.push_back(tpoint[i-1] + norm2(knot[i], knot[i-1]));
        }
        length = tpoint[n-1] + norm2(knot[n-1], knot[0]);

        // Compute D-splines
        std::vector<struct point> xt, yt;
        for (int i = 0; i < n; i++)
        {
            struct point xt0{tpoint[i], knot[i].x}, yt0{tpoint[i], knot[i].y};
            xt.push_back(xt0);
            yt.push_back(yt0);
        }
        struct point xtn{length, knot[0].x}, ytn{length, knot[0].y};
        xt.push_back(xtn); yt.push_back(ytn);

        Bspline B1(S32, xt, Nature), B2(S32, yt, Nature);
        Bx = B1; By = B2;
    }

    // Calculate x = x(t)
    double Getx(const double & t) const
    {
        return Bx.solve(t * length);
    }

    // Calculate y = y(t)
    double Gety(const double & t) const
    {
        return By.solve(t * length);
    }

    // Get point by t
    struct point Getpoint(const double & t) const
    {
        struct point out{Getx(t), Gety(t)};
        return out;
    }
};





# else
// Do nothing.
# endif